/** @file essl_impl_impl.h
 * @brief Implementation of essl_impl-enabled FFT engines
 *
 */
#ifndef QMCPLUSPLUS_FFT_ESSL_IMPL_H
#define QMCPLUSPLUS_FFT_ESSL_IMPL_H
#include <essl.h>
#include <iostream>
#include <vector>

//template<class T>
struct fft_essl_plan_type
{
  double* aux1;
  double* aux2;

  fft_essl_plan_type():aux1(0),aux2(0){}

  ~fft_essl_plan_type()
  {
    clean();
  }

  void clean()
  {
    if(aux1) delete [] aux1; aux1=0;
    if(aux2) delete [] aux2; aux2=0;
  }

  void resize(int n1, int n2)
  {
    aux1=new double[n1];
    aux2=new double[n2];
  }
};

template <typename T> struct fft_essl_base {};

/** base class for FFT engines using essl library
 */
template<>
struct fft_essl_base<double>
{
  static const int FFTESSL_FORWARD=1;
  static const int FFTESSL_BACKWARD=-1;

  typedef double real_type;
  typedef std::complex<double> complex_type;
  //typedef fft_essl_plan_type<double> fft_plan_type;
  typedef fft_essl_plan_type fft_plan_type;

  int dims_;
  int howmany_;
  int naux1;
  int naux2;
  int forward_offset;
  int backward_offset;
  fft_plan_type forward_plan;
  fft_plan_type backward_plan;

  ///default constructor
  fft_essl_base() { }

  ///virtual destructor to clean up the plans
  virtual ~fft_essl_base() 
  { 
    forward_plan.clean();
    backward_plan.clean();
  }

  /* complex-to-complex */
  void create_plans(int dims_in, int howmany_in, complex_type* in, complex_type* out)
  {
    forward_plan.clean();
    backward_plan.clean();

    dims_=dims_in;
    howmany_=howmany_in;
    forward_offset=dims_;
    backward_offset=dims_;

    naux1=(dims_>2048)?static_cast<int>(20000.+2.28*dims_)+1:20000;
    naux2=std::max(naux1,(2*dims_+256)*std::min(64,howmany_));

    forward_plan.resize(naux1,naux2);
    backward_plan.resize(naux1,naux2);

    dcft(1,in,1,dims_,out,1,dims_,dims_,howmany_,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1,forward_plan.aux2,naux2);

    dcft(1,out,1,dims_,in,1,dims_,dims_,howmany_,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
  }

  void create_plans(int dims_in, int howmany_in, real_type* in, real_type* out)
  {
    forward_plan.clean();
    backward_plan.clean();

    dims_=dims_in;
    howmany_=howmany_in;

    forward_offset=dims_+2;
    backward_offset=dims_/2+1;
    naux1=(dims_>4096)?static_cast<int>(20000.+1.64*dims_)+1:25000;
    naux2=naux1;
    forward_plan.resize(naux1,naux2);
    backward_plan.resize(naux1,naux2);

    drcft(1,in,forward_offset,out,backward_offset,dims_,howmany_,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);

    dcrft(1,out,backward_offset,in,forward_offset,dims_,howmany_,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
  }

  inline void execute_forward(complex_type* in, complex_type* out)
  {
    dcft(0,in,1,dims_,out,1,dims_,dims_,howmany_,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);
  }

  inline void execute_backward(complex_type* in, complex_type* out)
  {
    dcft(0,in,1,dims_,out,1,dims_,dims_,howmany_,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
  }

  inline void execute_forward(real_type* in, real_type* out)
  {
    drcft(0,in,forward_offset,out,backward_offset,dims_,howmany_,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);
  }

  inline void execute_backward(real_type* in, real_type* out)
  {
    dcrft(0,in,backward_offset,out,forward_offset,dims_,howmany_,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
  }
};

template<>
struct fft_essl_base<float>
{
  static const int FFTESSL_FORWARD=1;
  static const int FFTESSL_BACKWARD=-1;

  typedef float real_type;
  typedef std::complex<float> complex_type;
  //typedef fft_essl_plan_type<float> fft_plan_type;
  typedef fft_essl_plan_type fft_plan_type;

  int dims_;
  int howmany_;
  int naux1;
  int naux2;
  int forward_offset;
  int backward_offset;
  fft_plan_type forward_plan;
  fft_plan_type backward_plan;

  ///default constructor
  fft_essl_base() { }

  ///virtual destructor to clean up the plans
  virtual ~fft_essl_base() 
  { 
    forward_plan.clean();
    backward_plan.clean();
  }

  void create_plans(int dims_in, int howmany_in, complex_type* in, complex_type* out)
  {
    forward_plan.clean();
    backward_plan.clean();

    dims_=dims_in;
    howmany_=howmany_in;
    forward_offset=dims_;
    backward_offset=dims_;

    naux1=(dims_>8192)?static_cast<int>(20000.+1.14*dims_):20000;
    naux1=(naux1/4+1)*4;
    //naux2=std::max(naux1,(2*dims_+256)*std::min(64,howmany_));
    naux2=naux1;

    forward_plan.resize(naux1,naux2);
    backward_plan.resize(naux1,naux2);

    scft(1,in,1,dims_,out,1,dims_,dims_,howmany_,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1,forward_plan.aux2,naux2);

    scft(1,out,1,dims_,in,1,dims_,dims_,howmany_,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
  }

  void create_plans(int dims_in, int howmany_in, real_type* in, real_type* out)
  {
    forward_plan.clean();
    backward_plan.clean();

    dims_=dims_in;
    howmany_=howmany_in;

    forward_offset=dims_+2;
    backward_offset=dims_/2+1;
    naux1=(dims_>16384)?static_cast<int>(20000+0.82*dims_):25000;
    naux1=(naux1/4)*4;
    naux2=(dims_>16384)?static_cast<int>(20000+0.57*dims_):20000;
    naux2=(naux2/4)*4;


    forward_plan.resize(naux1,naux2);
    backward_plan.resize(naux1,naux2);
    srcft(1,in,forward_offset,out,backward_offset,dims_,howmany_,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1, naux1 ,forward_plan.aux2, naux2, 0, 0);

    scrft(1,out,backward_offset,in,forward_offset,dims_,howmany_,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1, naux1, backward_plan.aux2 ,naux2, 0, 0);
  }

  inline void execute_forward(complex_type* in, complex_type* out)
  {
    scft(0,in,1,dims_,out,1,dims_,dims_,howmany_,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1, naux1, forward_plan.aux2, naux2);
  }

  inline void execute_backward(complex_type* in, complex_type* out)
  {
    scft(0,in,1,dims_,out,1,dims_,dims_,howmany_,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1, naux1, backward_plan.aux2, naux2);
  }

  inline void execute_forward(real_type* in, real_type* out)
  {
    srcft(0,in,forward_offset,out,backward_offset,dims_,howmany_,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2,0,0);
  }

  inline void execute_backward(real_type* in, real_type* out)
  {
    scrft(0,in,backward_offset,out,forward_offset,dims_,howmany_,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1,backward_plan.aux2,naux2,0,0);
  }
};


/** specialization for single-fft 
 */
template<typename PREC, typename T>
struct FFT_DFT<PREC,T,false,FFTESSL_ENG> 
: public fft_base, public fft_essl_base<PREC>
{
  typedef T value_type;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,1,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      fft_base::resize(dims,1);
      fft_essl_base<PREC>::create_plans(dims,1,in,in);
    }
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_size) fft_essl_base<PREC>::execute_forward(offset,offset);
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_size) fft_essl_base<PREC>::execute_backward(offset,offset);
  }
};


/** specialization for multiple-fft with essl-FFT on complex
 */
template<typename PREC, typename T>
struct FFT_DFT<PREC,T,true,FFTESSL_ENG> 
: public fft_base, public fft_essl_base<PREC>
{
  typedef T value_type;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,howmany,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      fft_base::resize(dims,howmany);
      fft_essl_base<PREC>::create_plans(dims,howmany,in,in);
    }
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    fft_essl_base<PREC>::execute_forward(in,in);
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    fft_essl_base<PREC>::execute_backward(in,in);
  }
};

#endif
